{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "demonstrated-player",
   "metadata": {},
   "source": [
    "# Problem 14.08 Required Compressor Power for Isothermal and Adiabatic Compression of a Gas Mixture (CO2, O2) Using the Ideal Gas Law"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "different-transition",
   "metadata": {},
   "source": [
    "A stream of 1000 mol/hour CO2 and 1000 mol/hour O2 is compressed from 290 K and 1 bar to 5 bar. Calculate the compression power for both adiabatic compression, and isothermal compression. The compression is reversible (assumed) in each case - no efficiencies are necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "assigned-employee",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "This is a straightforward calculation. Using Thermo, working with complicated mixtures can be about as easy as pure components - if binary interaction parameters are zero. In this case, we try to load a parameter from a sample ChemSep database, but no values are available.\n",
    "\n",
    "The values in that database are just a sample - it is entirely the user's responsibility to provide the correct data to Thermo. If garbage is put in, garbage will come out!\n",
    "\n",
    "The problem says to use the ideal-gas law, so we can do that too and see how the answers compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "expressed-cartoon",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The PR kij matrix is [[0.0, 0], [0, 0.0]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.constants import hour\n",
    "T1 = 290\n",
    "P1 = 1e5\n",
    "P2 = 5e5\n",
    "flow = 2000/hour # mol/s\n",
    "\n",
    "from thermo import ChemicalConstantsPackage, PRMIX, IGMIX, FlashVL, CEOSLiquid, CEOSGas\n",
    "from thermo.interaction_parameters import IPDB\n",
    "\n",
    "constants, correlations = ChemicalConstantsPackage.from_IDs(['CO2', 'O2'])\n",
    "kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')\n",
    "print(f'The PR kij matrix is {kijs}')\n",
    "\n",
    "eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas,\n",
    "                 kijs=kijs)\n",
    "liquid = CEOSLiquid(PRMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)\n",
    "gas = CEOSGas(PRMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)\n",
    "flasher = FlashVL(constants, correlations, liquid=liquid, gas=gas)\n",
    "zs = [.5, .5]\n",
    "\n",
    "liquid = CEOSLiquid(IGMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)\n",
    "gas = CEOSGas(IGMIX, HeatCapacityGases=correlations.HeatCapacityGases, eos_kwargs=eos_kwargs)\n",
    "flasher_ideal = FlashVL(constants, correlations, liquid=liquid, gas=gas)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accompanied-encounter",
   "metadata": {},
   "source": [
    "### Adiabatic compression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "spatial-shadow",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shaft power with Peng-Robinson is 2632.7895 W\n",
      "The shaft power with ideal-gas is 2639.9248 W\n"
     ]
    }
   ],
   "source": [
    "# Solve with Peng-Robinson\n",
    "state_1 = flasher.flash(T=T1, P=P1, zs=zs)\n",
    "state_2 = flasher.flash(S=state_1.S(), P=P2, zs=zs)\n",
    "shaft_duty = (state_2.H() - state_1.H())*flow\n",
    "\n",
    "print(f'The shaft power with Peng-Robinson is {shaft_duty:.4f} W')\n",
    "\n",
    "\n",
    "state_1 = flasher_ideal.flash(T=T1, P=P1, zs=zs)\n",
    "state_2 = flasher_ideal.flash(S=state_1.S(), P=P2, zs=zs)\n",
    "shaft_duty = (state_2.H() - state_1.H())*flow\n",
    "print(f'The shaft power with ideal-gas is {shaft_duty:.4f} W')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enabling-process",
   "metadata": {},
   "source": [
    "### Isothermal Compression\n",
    "\n",
    "This problem is more interesting, because there is the cooling duty as well as the compressing duty.\n",
    "\n",
    "From theory, in an ideal gas, the cooling duty will be exactly equal to the compressing duty.\n",
    "\n",
    "For a real-gas, it will be different as enthalpy is pressure-dependent.\n",
    "\n",
    "In both cases, the evaluation of the following integral is required.\n",
    "\n",
    "$$ \\text{duty} = \\text{flow} \\int_{P1}^{P2} V \\partial P $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "everyday-flesh",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shaft power with ideal-gas is 2155.9263 W\n",
      "The cooling duty with ideal-gas is 2155.9263 W\n",
      "The shaft power with Peng-Robinson is 2139.44610002 W\n",
      "The cooling duty with Peng-Robinson is 2192.57596810 W\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "def to_int(P, flasher):\n",
    "    state = flasher.flash(T=T1, P=P, zs=zs)\n",
    "    return state.V()\n",
    "shaft_duty = cooling_duty = quad(to_int, P1, P2, args=(flasher_ideal,))[0]*flow\n",
    "\n",
    "print(f'The shaft power with ideal-gas is {shaft_duty:.4f} W')\n",
    "print(f'The cooling duty with ideal-gas is {cooling_duty:.4f} W')\n",
    "\n",
    "entry = flasher.flash(T=T1, P=P1, zs=zs)\n",
    "exit = flasher.flash(T=T1, P=P2, zs=zs)\n",
    "\n",
    "shaft_duty = quad(to_int, P1, P2, args=(flasher,))[0]*flow\n",
    "cooling_duty = shaft_duty - (exit.H() - entry.H())*flow\n",
    "\n",
    "print(f'The shaft power with Peng-Robinson is {shaft_duty:.8f} W')\n",
    "print(f'The cooling duty with Peng-Robinson is {cooling_duty:.8f} W')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "broadband-excerpt",
   "metadata": {},
   "source": [
    "The above shows the numerical integral calculation. That is the correct formulation.\n",
    "\n",
    "However, it can be a little unintuitive. We can contrast this with another calculation - a series of tiny isentropic compression, then cooling steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "noticed-shaft",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shaft power is 2322.61227046 W\n",
      "The cooling duty is 2375.74213854 W\n"
     ]
    }
   ],
   "source": [
    "cooling_duty = 0\n",
    "compressing_duty = 0\n",
    "increments = 3 # Number of increments\n",
    "dP = (P2 - P1)/increments\n",
    "old_state = entry\n",
    "for i in range(increments):\n",
    "    P = P1+(i+1)*dP\n",
    "    \n",
    "    # Compress another increment of pressure\n",
    "    new_compressed_state = flasher.flash(S=old_state.S(), P=P, zs=zs)\n",
    "    compressing_duty += (new_compressed_state.H() - old_state.H())*flow\n",
    "    \n",
    "    # Cool back to T1 at new pressure\n",
    "    new_cooled_state = flasher.flash(T=T1, P=P, zs=zs)\n",
    "    cooling_duty += (new_compressed_state.H() - new_cooled_state.H())*flow\n",
    "    \n",
    "    old_state = new_cooled_state\n",
    "\n",
    "print(f'The shaft power is {compressing_duty:.8f} W')\n",
    "print(f'The cooling duty is {cooling_duty:.8f} W')"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}